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Abstract. A family of triaxial Stackel potential-density 
pairs is introduced. With the help of a Quadratic Pro- 
gramming method, a linear combination of potential- 
density pairs of this family which fits a given projected 
density distribution can be built. This deprojection strat- 
egy can be used to model the potentials of triaxial el- 
liptical galaxies with or without dark halos. Besides, we 
show that the expressions for the Stackel triaxial density 
and potential are considerably simplified when expressed 
in terms of divided differences, which are convenient nu- 
merically. We present an example of triaxial deprojection 
for the galaxy NGC 5128 whose photometry follows the 
de Vaucouleurs law. 
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1. Introduction 

As soon as one realized that elliptical galaxies were not 
spheroids flattened by rotation, a growing list of unex- 
pected photometric and kinematical observations has led 
astronomers to develop a more elaborate picture of these 
objects. While triaxiality provides a plausible explanation 
of these new findings, it also raises a number of questions 
such as the determination of the intrinsic shapes of these 
stellar systems. 

The deprojection problem consists in finding the in- 
trinsic mass (or light) density from its projection onto the 
plane of the sky. It is well known that the knowledge of the 
surface brightness distribution does not define uniquely 
the intrinsic three-dimensional light distribution and the 
orientation of the galaxy. Even in the axisymmetric case 
the deprojection is degenerate: one can construct many 
density distributions that project to the same photomet- 
ric distribution (see e.g. Gerhard & Binney 1995). Features 



such as the presence of lanes or disks in some galaxies have 
been exploited to further constrain the problem. Recently 
Statler (1994) proposed a new method for constraining the 
intrinsic shapes of elliptical galaxies using information on 
their apparent shapes and velocity fields. 

In this paper, we shall describe a technique that de- 
termines an intrinsic mass density distribution of a galaxy 
consistent with a given projected mass density. Similarly, 
the method allows to determine the intrinsic light distri- 
bution of a galaxy given its surface brightness distribution. 
This is the first step towards a triaxial dynamical model. 
Stackel potentials are particularly suited to this kind of 
modeling since the investigation of the dynamics is es- 
sentially analytical. A Stackel triaxial mass model can be 
constructed by specifying a density profile along the short 
axis and an ellipsoidal coordinate system. Simple examples 
of such models are given in de Zeeuw, Peletier & Franx 
(1986). However, the determination of the potential gen- 
erally requires the evaluation of one-dimensional quadra- 
tures. Another approach consists in specifying ellipsoidal 
coordinates and a simple form for the Stackel potential. 
The mass density can be calculated by means of the gener- 
alized Kuzmin's formula (de Zeeuw 1985b) which involves 
straightforward derivatives of a one-dimensional function. 
We shall use this method in order to have a simple ana- 
lytical expression for the potential. 

The basic properties of Stackel models are studied in 
Sect. 2. Section 3 presents a family of Stackel potential- 
density pairs that can be used as building blocks for the 
construction of triaxial mass models. Examples of pro- 
jected and spatial density distributions are shown and 
then we describe the deprojection method. In Sect. 4, 
we present a fit to a triaxial modified Hubble model us- 
ing these potential-density pairs. Section 5 deals with the 
application of the deprojection method to the elliptical 
galaxy NGC 5128 (Centaurus A) whose photometry fol- 
lows the de Vaucouleurs law. Our conclusions are given in 
Sect. 6. 
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2. Triaxial Stackel models 

A sequence of potential-density pairs can be generated 
by repeated application of an operator upon a given 
potential-density pair. Such a method has been used by 
de Zeeuw and Pfenniger (1988, hereafter referred to as 
ZP88) to build a class of triaxial potential-density pairs 
(V n ,Pn)i amongst which the family (V^,/^) is of Stackel 
form. In the next subsection, we present some properties of 
the Stackel models (for a detailed discussion see de Zeeuw 
1985a, ZP88). 

2.1. General properties 

Let (x, y, z) be cartesian coordinates and (A, p, v) be ellip- 
soidal coordinates defined according to de Zeeuw (1985a). 
In these coordinates, the Stackel potential, for which the 
equations of motion separate, has the general form 



V(X, p, v) = g x F{\) + 9fl F(p) + g v F{y), 
where F is an arbitrary function and 
(X + a)(X + [3) 



fix 



{eye). 



{X-p)(X-v) 
a, (3 and 7 are negative constants and 
-7 < v < -(3 < p < ~a < A. 



(1) 



(2) 



(3) 



Similar expressions for g M and g v are found by a cyclic 
permutation A — * p — > v — > A. Furthermore, we have 
< g\, sy, , g v < 1 and 



5a + g^ + g v = 1. 



(4) 



The triaxial density follows from Poisson's equation and 
can be written (see ZP88) as 

p(X, p, v) = g\V{X) + g 2 ^'{p) + gl^\v) + (5) 

2g x g f2 W[X, p] + 2 9tl g^[p, v] + 2g v g x ^[u, A], 

where ^[ti^] denotes the first order divided difference 
of the function ^(r) i.e. 



*[ri,T 2 ] 



fr(Tl) - #(T 2 ) 

n - ti 



(G) 



and ^'(t) is the derivative \V[t, t]. The relation that con- 
nects *(t) and F(t) is (ZP88): 



2ttG*(t) = 2(t + j)F'(t) - F(t) + 

( j ( F(r)-F(-a) | F(r)-F(-/3) 
V t + a t + (3 



(7) 



Given a function F(t), one can calculate the Stackel po- 
tential with Eq. (|TJ), and the density using Eq. (§). The 
density profile along the z-axis is the function ^'(r) with 
r = z 2 — 7. The triaxial Stackel density is fully determined 
by the specification of the density on the short z-axis; it 
can be written as a weighted sum of the density at six 
particular points of the z-axis (de Zeeuw 1985b). 



2.2. Expressions for the Stackel potential and density in 
terms of divided differences 



The use of divided differences considerably simplifies the 
formulas of the Stackel triaxial density and potential. The 
divided difference of order (n — 1) of the function U(t) is 
a function of divided differences of order {n — 2) and is 
given by 



U[ti,t 2 , ■ 



U[n,T 3 , 



U[t 2 ,t 3 , 



Tl - T 2 



(8) 



The second order divided difference of a general function 
U(t) can be written explicitly as 



U[X,p, v] 



U{X) 



{X-p){X-v) 



{eye). 



(9) 



With Eqs. (g) and (|]J), we can express the potential V as 
the second order divided difference 



V(X,p,v) = U[X,p,v] 
with 

U(t) = (t + a)(r + @)F{t). 



(10) 



(11) 



The density p(X, p, v) can be written in terms of di- 
vided differences using either the function \&(t) whose 
derivative is the density profile along the z-axis, or us- 
ing the basis function F(t) that defines the potential V. 
The detailed proofs are given in Appendix B and Ap- 
pendix C respectively. To establish these results, we shall 
use the potential-density pair (pi, V\) introduced in ZP88. 
Its main properties are given in Appendix A. In this sub- 
section, we only present the relevant results. 

The triaxial Stackel density p(X, p, v) can be expressed 
as the 5 th -order divided difference of the function H{t) 



p(X, p, v) = H[X, p, v, A, p, v\ 
with 

H{t) = (r + a) 2 (r + /3) 2 *(r). 



(12) 



(13) 



It can also be written as the 5 th -order divided differ- 
ence of the function i?(r) 



p(X,p,v) = R[X, p,, v, A, p, v\ 
with 



■KGR{r) = \a{r)\-^ 



t + 7 



F(r) 



where 

a{r) = (r + a)(r + 0)(t + 7). 



(14) 



(15) 



(16) 
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3. The deprojection strategy 

3.1. Basis component 

A Stackel triaxial mass model is completely determined by 
the choice of ellipsoidal coordinates and a function F(t). 
So as to have simple analytical expressions for the poten- 
tial and the density, we choose F(r) to be an elementary 
function of r. 

The spherical Henon's isochrone model (1959) is of- 
ten used in the construction of models of stellar systems. 
More realistic models that account for the flattening of 
the potential can be produced with axisymmetric general- 
izations (see Dejonghe & de Zeeuw 1988; Evans, de Zeeuw 
& Lynden-Bell 1990). Axisymmetric Stackel models that 
reduce to the isochrone in the spherical limit can be con- 
structed using 



In this section, we present three examples of spatial 
and projected density distributions for a normal ellipti- 
cal (Fig. |l|), a discy (Fig. ||) and a boxy (Fig. ||) mod- 
els. The units are arbitrary. The intrinsic long, interme- 
diate and short axes are denoted x, y and z respectively. 
The models have (-a , -(3 , -7) = (7.3 , 3.6 , 0.8). The pa- 
rameters (d , p , s) of the components are (5 , 0.7 , 0.7) for 
the elliptical model, (4,0.5,1) for the discy model and 
(5 , 0.93 , 0.44) for the boxy model. 



F(r) 



GM 



(17) 



The associated potential is the thoroughly studied 
Kuzmin-Kutuzov potential (see e.g. Kuzmin 1956, 
Kuzmin & Kutuzov 1962, Dejonghe & de Zeeuw 1988, Bat- 
sleer & Dejonghe 1993). Using the above function F(r), a 
triaxial generalization can also be calculated with Eqs. (El) 
and (||) or, alternatively, Eqs. (]ll) and (|l5|). Expressions 
for the density and the potential are given in ZP88. 

We choose a three-parameter basis function of the form 



F(r) = 



GM 

(d+rpy 



(18) 



where d, p and s are real parameters. The triaxial 
isochrone has p = 0.5, s = 1 and d = yj—a.. 

3.2. Some examples 

Kinematic studies of elliptical galaxies have revealed a 
very diverse and complicated nature of these objects. 
There has been marked improvement in the photometric 
observations with the advent of CCDs and it appears that 
ellipticals also span a wide range of photometric proper- 
ties. The isophotal shape is mainly elliptical, but depar- 
tures from perfect ellipses are often detected, as well as 
isophote twists (see e.g. Bender, Dobereiner & Mollenhoff 



The projected mass densities of the basis components 
of the form (|lj) have elliptical surface isodensities that 
can exhibit deviations from pure ellipses such as boxiness 
and diskiness. Thus these components may be relevant 
to approximate the photometry of elliptical galaxies. It 
is well known that Stackel potentials cannot produce pro- 
jected densities with isophote twists (Franx 1988), but this 
is often a second order effect compared to the ellipticity 
variation. 



Fig. 1. Contour map of the spatial mass density in the planes 
(top left) x = 0, (top right) y — and (bottom left) z — for 
an elliptical model. Bottom right: contour map of the projected 
mass density. The contour step is 0.65 magnitude 



3.3. The method 

We use a Quadratic Programming (QP) method (De- 
jonghe 1989) to fit to a given photometric data set a linear 
combination of basis functions with the constraint that 
the density must be positive. Given a function F(t) of 
the form (|l8|), the calculation of the Stackel potential and 
density is straightforward as it involves only the evalua- 
tion of elementary functions. These basis components can 
produce physical density distributions (i.e. positive every- 
where in real space) as well as unphysical ones, depend- 
ing on the choice of the ellipsoidal coordinates and the 
parameters d, p and s. Generally, the projected density 
distribution has to be calculated numerically. 

As showed by Kuzmin (1956) for axisymmetric mass 
models and subsequently generalized for triaxial models 
by de Zeeuw (1985b), the so-called Kuzmin's theorem 
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states that a triaxial mass density built with a Stackel 
gravitational potential is everywhere positive if the density 
on the z-axis is positive. Therefore, the positivity check of 
the triaxial density p{\, p, v) reduces to a one-dimensional 
check on the density along the z-axis. 



4. Triaxial modified Hubble model 

The de Vaucouleurs and the Hubble laws are among the 
most commonly used fitting profiles for the photometry of 
ellipticals. The density radial profile of a modified Hub- 



ble model p ~ (1 



" 3//2 decreases proportional to 



at large radii. A separable triaxial generalization of the 
modified Hubble model has been studied by de Zeeuw, 
Peletier & Franx (1986). The density along the z-axis for 
a modified Hubble model with core radius c can be written 



Fig. 2. Contour map of the spatial mass density in the planes 
(top left) x — 0, (top right) y = and (bottom left) z — for a 
discy model. Bottom right: contour map of the projected mass 
density. The contour step is 0.6 magnitude 



Poc" 



(z 2 + C 2 ) 3 / 2 



By choosing 7 
it becomes 



(19) 



and using z 2 = t + 7 on the z-axis, 



*'(tHpo( — 



3/2 



The function W(t) is the primitive 



*(r) = j ^'(a)da = 2p o7 ( 



(20) 



(21) 



The triaxial density can be calculated with Eqs. (|T^) and 

Using basis components of the form (tm) , we can fit the 
density on the z-axis ^'(r) of a modified Hubble model. 
As the triaxial density is completely determined by the 
ellipsoidal coordinate system and the function \P'(t), we 
can therefore produce a fit to a triaxial Hubble model. 
A wide range of triaxiality is allowed through the choice 
of the ellipsoidal coordinates. In Fig. |[ we present the 
spatial and projected mass densities of our best fitting 
model with (— a,— (3,— 7) = (4,2.6,1). The model is 
fully triaxial, with a triaxiality parameter T defined as 
T = (A 2 — B 2 )/(A 2 - C 2 ) (with A, B and C the long, in- 
termediate and short axis lengths of the density) of ~ 0.5. 
The relative difference between the function \&'(t) and the 
fit is smaller than 1 % out to z ~ 40 c. 



Fig. 3. Contour map of the spatial mass density in the planes 
(top left) x — 0, (top right) y = and (bottom left) z — for a 
boxy model. Bottom right: contour map of the projected mass 
density. The contour step is 0.6 magnitude 



5. Deprojection of a de Vaucouleurs photometry 

Centaurus A (NGC 5128) is a giant elliptical galaxy with 
a conspicuous dust lane lying along its photometric mi- 
nor axis. A photometric study of Cen A by Dufour et al. 
(1979) from photographic plates showed that the light dis- 
tribution of Cen A follows the de Vaucouleurs law at radial 
distances from 2 arcmin to 8 arcmin. Furthermore, they 
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Fig. 4. Contour map of the spatial mass density in the planes 
(top left) x — 0, (top right) y — and (bottom left) z — 
for the best fitting model to a triaxial modified Hubble model. 
Bottom right: contour map of the projected mass density (solid 
line) of the fit. The contour step is 0.5 magnitude. The dashed 
line is the projected mass density of the modified Hubble 
model. The direction of projection is the j/-axis 



found that the V surface brightness distribution at radii 
from 4 arcmin up to 8 arcmin is consistent with Cen A 
being an E2 galaxy. The isophotes are elliptical in shape 
and are quite round at the center but become more flat- 
tened with increasing radius. The ellipticity defined as 
e = 1 — b/a (with a and b the apparent major and minor 
axis lengths) increases from e = 0.07 at r = 2.6 arcmin 
to e = 0.26 at r = 9 arcmin. Following Hui et al. (1993), 
we adopt a distance of 3.5 Mpc to Cen A, so that 1 ar- 
cmin corresponds to 1.02 kpc. According to Dufour et al. 
(1979), the effective radius of the de Vaucouleurs law is 
r e = 5.18 kpc. 

We build a model for the photometry of Cen A that 
follows the de Vaucouleurs law and reproduces the ob- 
served flattening of the light distribution out to 8 arcmin. 
Beyond this radius, we choose a constant value of the el- 
lipticity e = 0.26 and a wide range of mass density profiles 
is allowed that can account for a possible increase of the 
mass-to-light ratio with radius. We adopt the observer's 
viewing direction determined by Hui et al. (1995). 

We use the QP method described in Sect. 3 to fit a 
Stackel model to the de Vaucouleurs photometry of Cen A. 
Given ellipsoidal coordinates, we find that models ranging 
from quite oblate to fully triaxial can fit the data. In gen- 
eral, our QP models consist of less than ~ 10 components. 



Fig. 5. Contour map of the spatial light density in the planes 
(top left) x — 0, (top right) y = and (bottom left) z = for 
the best fit model to Cen A's photometry. Bottom right: con- 
tour map of the projected light density (solid line). The dashed 
line is the de Vaucouleurs photometry model. The contour step 
is 0.5 magnitude 



Adding new components does not improve the fit signifi- 
cantly. We present a model with a triaxiality parameter T 
of ~ 0.5. Contours of the spatial density in the principal 
planes are shown in Fig. |[ The contours of the projected 
density (solid line) of the QP best fit are compared to the 
contours of the de Vaucouleurs photometry model (dashed 
line) . The residual differences of 92 % of the pixels inside 
the region limited by the isophote of major axis length of 
35 arcmin are smaller than 0.1 magnitude. Figure |^ shows 
the profiles of the QP model and the de Vaucouleurs pho- 
tometry model along the photometric major and minor 
axes and the apparent axis ratio as a function of radius 
for the two models. The components allow to reproduce 
the de Vaucouleurs law in the range ~ 0.5 r e to ~ 6r e . 
Taking into account the fact that the mass-to-light ratio 
is likely to vary a bit with radius in this galaxy, there is 
no particular need to attempt a better fit if a deprojected 
mass model and hence potential are desired. The fit yields 
a sufficient approximation of the galaxy's potential. In the 
inner region (r < 0.5 r e ), the fit is less steep than the de 
Vaucouleurs law (see Fig. ^). Dynamical models produced 
with this potential may depend on three integrals, which 
will be only approximative in the very center. 
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exhibits unambiguous signatures of triaxiality. Using a 
Quadratic Programming method, we find a linear com- 
bination of density distributions that fits a model of the 
surface brightness of this E2 galaxy with a total density 
which is positive everywhere. The de Vaucoulcurs photom- 
etry model of Cen A is well reproduced with these com- 
ponents. Mass-to- light ratio variations could be included 
in the projected distribution and the same method would 
then produce spatial mass models with dark matter. It 
suggests that Stackel mass models may be relevant to the 
description of galactic potentials which represents a first 
step towards Stackel triaxial dynamical models of stel- 
lar systems. Abel components propounded by Dejonghe & 
Laurent (1991) may prove powerful in this context; work 
along this line was first carried out by Dejonghe (1992), 
subsequently by Zcilinger ct al. (1993) and further inves- 
tigation is in progress. 

Acknowledgements. We thank the referee E. Emsellem for his 
critical comments which helped us improve the presentation. 
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Fig. 6. Top: profiles of the projected density of our QP best 
fit model to Cen A's photometry along the photometric major 
(bottom solid line) and minor (top solid line) axes. The dashed 
lines are the corresponding profiles of the de Vaucouleurs pho- 
tometry model. Bottom: apparent axis ratio as a function of 
radius for the QP model (solid line) and for the de Vaucouleurs 
photometry model (dashed line) 



6. Conclusions 

A family of triaxial Stackel potential-density pairs is pre- 
sented. It includes as a special case the triaxial generaliza- 
tion of Henon's isochrone in ellipsoidal coordinates. This 
family allows the construction of triaxial mass (or light) 
models of galaxies with or without dark halos. A large 
variety of intrinsic shapes is provided by the choice of el- 
lipsoidal coordinates and by the components. This diver- 
sity is also reflected in the projected densities which show 
elliptical, box-like or disc-like isophotes with ellipticities 
changing as a function of radius. 

These potential-density pairs can be used as building 
blocks for realistic Stackel models of triaxial potentials 
in elliptical galaxies. This is first tested with the ellip- 
tical galaxy Centaurus A (NGC 5128) whose kinematics 



A. The potential-density pair (p 1 ,V 1 ) 

Let the potential V\ be of the general form 

V 1 (\,»,v)=F 1 (\)+F 1 (li)+F 1 (v). (Al) 

The Laplace operator V 2 in ellipsoidal coordinates is (see 
e.g. de Zeeuw 1985b, Eq. (20)) 

V 2 =V!+V 2 +V 2 (A2) 



with 



(A-mKA-,,) { 2a ^ + a '^X 



where 

a(r) = (r + a)(T + /3)(r + 7). 



(eye), (A3) 



(A4) 



The density p\ that follows from Poisson's equation can 
be written as (Eq. (3.9) in ZP88) 



Pl (\,p,v) = 5A *i(A) + (eye). 



(A5) 



The relation between *i(r) and F\(t) is written explicitly 
in ZP88 (Eq. (3.7)), as well as the relations that connect 
the Stackel pair (p, V) and the pair (p\, Vi). In particular, 
we have 



p={A + l)p u 
F(t) = (t + 7 )F{(t), 

*(t) = (t + 7)*i(t) 
with 

d_ 

l 9A 



A = (A + j)gx— + (p + 7)5^ + (" + 7)3- — 



d_ 

dv 



(A6) 
(A7) 
(A8) 



(A9) 



and X is the identity operator. 
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B. The triaxial density as a function of \P(t) 

By noting that H [n, n, t 2 , t 3 ] = ^H[ti, t 2 , t 3 ], we have 



d d d 

H[X,p,v, A,/i, i/] = — — — ff[A,pt,i/]. 

c/A Of 

Let if(r) be the function 
H{t) = (T + a) 2 (r + /3) 2 *(r). 



Then Eq. (Bl) can be written 
H[\,n,v, A,/i, v\ = 



(Bl) 

(B2) 
(B3) 



d_ 

dX 



(A + a )(A + /?)*(A) — — g x 
Op 01/ 



or, with Eq. (|), 



(eye), 



(B4) 



iJ[A,/i,^A,/i,z/]= fl ^'(A)+2 flA *(A)^ + (cyc). (B5) 



9A 



Using Eq. (|]), we have 

dflA = 9^ 9_u_ 

dX (A - fi) (A - i 



(eye). 



(B6) 



Substitution of Eq. (B6) in Eq. (B5) and comparison with 
Eq. (||) prove that the density p can be expressed as the 
5 th -order divided difference 



P(A, p, v) = H[X, p, v, A, p, u] 



(B7) 



with H{t) given in Eq. ( |B2|) . By specifying a function 
^(t) or the density profile on the z-axis ^'(t), one can 
easily calculate the density with Eq. ( B7 ) . One can notice 
that any 4 th -order polynomial added to the function H(t) 
yields the same 5 th -order divided difference as H(t). 

C. The triaxial density as a function of F(t) 

One can also calculate the density from the potential by 
specifying the function F(t) instead of ^(r). We show 
that the triaxial density can be written as a 5 th -order di- 
vided difference of a function R(t) that depends only on 
F(r). 

By direct application of Poisson's equation, using 
Eqs. (Al), ( |A3| ) and (^|), we can write the density p\ as 



pi(X,p,v) = i?i[A, p, v] 
with 

kGR x {t) = a{r)F^r)+ l -a'{T)F[{T) 



a(r) 



dr 



[VWY\F[(r)}. 



(CI) 

(C2) 
(C3) 



Now we prove that the density p can be written as the 
5 th -order divided difference 



p(X,p,v) = R[X, fi, v, A, fi, v] 
where 

R{t)=cl(t)R 1 {t). 



(C4) 



(C5) 



Using Eq. (Bl), we can write 
R[X,p, v,X,ii,v] 



d d d 
dX dp dv 



d_ 

dX 



a(X)Ri(X) 
a{X)R 1 (X) 



{X-py{X-vf 



Since 



a(A) 



{X-p)(X-v) 
we obtain 



(A + 7).9a (eye), 



R[X,p,,v,X,p,,v] = {X + -y)g\ 



dX 



Ri(X) 



(X-p)(X-v) 



(eye) 
(eye). (C6) 

(C7) 

(C8) 



Ri(X) 



-[{X + j)gx] + {cyc), 



(X-p){X-is) dX 1 
With Eqs. (|) and (^6|) , one can establish the identity 



^[(A + 7).9a] = 1 + ^.9, + ^.9, (eye), (C9) 



which, upon substitution of Eq. (C9) in Eq. (C8), yields 

R[X,p,,v,X,p,,v] = (A + l)R 1 [X 1 p,v\ = (A + l)pi, (CIO) 

identical to Eq. (|A6|), q.e.d. 

Thus, with Eq. (|A7|), the density p can be written as 
the 5 th -order divided difference (Eq. (p4)) with 



hGR{t) = |a(r)|3 



dr 



T + 7 



F(r) 



(Gil) 



and a(r) is given in Eq. (A4). The function R(t) differs 
from the function H (r) only by a polynomial function of 
the 4 th -order in r. 
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